An elementary mode coupling theory of random heteropolymer 

dynamics 

(glass transition / protein folding / Langevin equation ) 
Classification: Physics, Biophysics 

Shoji Takada, John J. Portman, and Peter G. Wolynes 
School of Chemical Sciences and Department of Physics, University of Illinois, Urbana, IL, 61801 

Contributed by Peter G. Wolynes 
Abbreviation: MCT, mode coupling theory ; RHP, random heteropolymer ; 
DFI, dynamic functional integral. 
(February 1, 2008) 



Abstract 

The Langevin dynamics of a random heteropolymer and its dynamic glass 
transition are studied using elementary mode coupling theory. Contrary to 
recent reports using a similar framework, a discontinuous ergodic-nonergodic 
phase transition is predicted for all Rouse modes at a finite temperature Ta- 
For sufficiently long chains, Ta is almost independent of chain length and is in 
good agreement with the value previously estimated by a static replica theory. 
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The self-organization of evolved biopolymers such as foldable proteins ultimately depends 
upon the chain dynamics of heteropolymers. The modern statistical mechanical view using 
energy landscape ideas focuses on the analogy of folding with phase transitions in finite size 
systems and exploits to a large extent our understanding of the quasi-thermodynamic and 
static features of both regular systems and frustrated disordered systems such as spin glasses 
IJ. The connection with detailed analytical theories of the time dependence of fluctuational 
motions of a randomly interacting chain has only recently received attention ^] . 

For completely random heteropolymers (RHP) the quasi-thermodynamic analyses ||[7j 
using replica techniques give results parallel to those for other systems with random first 
order phase transitions such as the mean field Potts spin glasses [|J. An important feature of 
theories of random first order transitions is the presence of two transitions; one is static, while 
the other is dynamic and generally occurs at a higher temperature. The dynamical transition 
signals a crossover to a motional mechanism involving activated motions ||,[| (We denote 
this dynamical transition temperature T\). For RHP's two different kinds of approximate 
quasi-static theories based on replicas ||[7| and on the generalized random energy model 
]iT| that takes into account correlation in the energy landscape yield these two transitions. 
Replica based techniques also yield estimates for the free energy barrier heights of activated 
motions between the two characteristic temperatures f||7]]. Purely dynamical theories based 
on mode coupling ideas generally yield a transition in harmony with the quasi-static analysis 
of the dynamical transition which is in some sense a spinodal. While mode coupling theory 
(MCT) is perturbatively well defined for spin systems with long range random interactions 
|TT|,|9[| , there are various versions developed in the theory of fluids [0 and polymers |H| that 



are forced to make uncontrolled approximations. On the basis of two such calculations, the 
validity of the emerging picture of the dynamics and the energy landscape of RHP's have 
been questioned J4|. 

Roan and Shakhnovich derive a mode coupling theory of the RHP and explicitly solve 
it for a polymer in a good solvent [Q. It is no surprise that an uncollapsed chain has no 
dynamical transition, but the authors further claim that this is actually a structural fea- 



ture of the RHP mode coupling equations independent of the state of collapse. Thirumalai 
et al. || derive another set of mode coupling equations for a somewhat different model 
that exhibits no static replica symmetry breaking |TJ| and conclude there is a dynamical 



transition with a transition temperature that depends on the length scale of the motional 
mode considered ||. A numerical treatment of self-consistent equations for heteropolymer 
collapse, on the other hand, does show evidence for a dynamic freezing transition ||. The 
inconsistency of these results with each other and in the first two cases with replica calcula- 
tions is disturbing. The technical intricacy of these mode coupling calculations is a barrier 
to understanding their inconsistencies. We have therefore derived mode coupling equations 
using 'elementary' methods like those used decades ago in the theory of critical phenomena 
|12|| . The resulting equations differ in some respects from those of the earlier workers but 



clearly yield transitions in harmony with the quasi-static results. These equations also lead 
to an understanding of how the dynamical freezing depends on chain length and state of 
collapse as well as how the dynamics varies for different modes of chain motion. 
We consider a standard Hamiltonian for iV interacting beads 

JV 

H = k B Tj2(r l+1 - r,) 2 + 1/2 ^^(Ary) + - J2Mt) ■ r ^)> (!) 

i=l i^j 

where Tj are the bead locations, by is chosen as Gaussian random with mean bo and variance 
b 2 , Ar-jj = r.j — rj, and u(r) is the two-body interaction chosen as a Gaussian exp — r 2 /cr 2 . V e ^ 
is the excluded volume term which usually contains three body terms, but it is enough, at 
the level of the present analysis for phantom chain, to replace it with an effective Gaussian 
confinement term V con f = kBTBJ2 r h with the constant B chosen so that the radius of 
gyration R g becomes the physically required value determined by the packing fraction, t] 
|15|j . hj is the external force introduced for convenience in the derivation of the response 



function. For technical simplicity, we adopt a ring polymer model where rjy+i = r^.. The 
Langevin equation for the beads is 

r- 1 ^ = -d(3H/dn + & (2) 
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Here is a Gaussian random force for which the first two moments are given by = 

and (£i(t)£j(t')) = 2T~ 1 5ij5(t—t')l (1 is the 3 by 3 unit matrix), T determines the microscopic 
time scale (it will be set to unity) and (3 = 1/ksT as usual. Since our main interest is in 
collapsed states where hydrodynamic effects are less important, we ignore them for simplicity 
although their inclusion is not very hard. A more essential simplification of our treatment 
is the neglect of entanglement effects. An uncrossable chain will possess more friction than 



the phantom chain treated here |L3|], modifying the dynamical transition temperature. 

A standard tool for dealing with the average over the quenched randomness in the 
Langevin dynamics has been the dynamic functional integral (DFI) formalism fTT|,|9|j. Al- 
though this approach is generally accepted for the infinite range spin model, direct applica- 
tion of the formalism to RHP models is questionable for two reasons: l)The heteropolymer 
interaction is short range so that the mean field approximation is not exact; in such a situa- 
tion, there exists some ambiguity with the DFI formalism in its use of steepest descents. 2) In 
contrast to spin models, the Jacobian appearing in the DFI formalism of the heteropoly- 
mer model depends on the randomness. This point was neglected in previous treatments. 
The usual formal identity for averaging over randomness in the Hamiltonian for a spin-glass 
model is no longer precisely true. Thus, we take an alternative and simpler route here; we 
use a perturbation theory in the randomness bij along with a self-consistent prescription that 
corresponds to a particular resummation of higher order terms. In the case of the infinite 
range p-spin model (including the Sherrington-Kirkpatrick model) , this procedure still leads 
to exactly the same results as the DFI formalism H|. 



We sketch the perturbation scheme briefly since details will be published elsewhere |L7 . 
The Langevin equation, Eq (|2|), can be expressed as 

r,(t) = £ G ,u * (& + hi - (3/2^- M ( Ar *;)) > ( 3 ) 

i 1 kj 

where * represents the time convolution and Go^j(t — t') is the 0-th order response function 
d(r\°\t)} /dhj(t'). Inserting the 0-th order solution = J2 Co,/i * £i into the right hand 
side of eq. (^) gives the first order solution = — f3/2J2i * 0/ dr[ ^ Y^kj bkjU ( Ar i°' ) ) • 



We insert this once again into eq.(^|), yielding the second order solution. After taking an 
average over 5^, we multiply the resulting equation by Gq 1 which yields 

r-%ri 2) = - £ K&f (t) +/i]T dt'AG 0)ij (t - Wu (Arg>(t)) ■ V« (Arg>(f )) + 6(f), 

(4) 

where l£y is the harmonic constant matrix including the confinement term, 
i^l = d 2 /dr l dr J k B T[E(r i+1 -r t ) 2 + BEr 2 }, fi = ((3b) 2 , and AG 0iij = G 0>ii - G 0>ij - 
Go,ji + ^oj'j- The random noise correlation can also be calculated by perturbation theory. 
With the use of the relation G\ m = Gu * (4>i4>j) * Gmj, we can derive the expression for 
the correlations up to second order, C^ 2 \ which includes {<j>iif)<j>j{lf))^ ■ 

A part of higher order terms can be taken into account by employing a self-consistent 
prescription. To this end, we first expand the memory kernel of eq.((4|) in Ar^t'). Then, 
in the spirit of Kawasaki's derivation of MCT of critical dynamics, we replace r- '* in the 
right hand side by i", and Go by the perturbed response function, G, giving a renormalized 
Langevin equation, 

r-Mi-i/di = - £ KijT^t) + /i E T dt'AG^t - t)dj(t - t')Ar l3 (t') + &(t), (5) 

• J —oo 
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where = (d [VVtt (AiY,-(t)) • Vw (Ar^ (£'))] /dAvijit')). In deriving this equation we used 
the isotropic symmetry of the model. The random force <pi(t) is also renormalized in the 
same way. In the expression for ((f)i(t)(f)j(t'))^ , replacement of by leads to the colored 
noise correlation function, 



(t)0j(O) = 2T- 1 15 ij 5{t-t')+ fi 



8 i:j Y,M ik {t-t')-M i:j {t-t') 



(6) 



where Mij = (Vu (Ar„(f)) Vw (Ar^ (£'))) is the force-force correlation on different beads. 
In calculating C and M we approximate the stochastic process as Gaussian. The ex- 
plicit expression for M becomes A%(t) = [(1 + sAC^O)) 2 - (sAC^(t)) 2 ] -5/2 AC^t), where 
ACij = Cu — Cij —Cji + Cjj and s = Ija 2 . For the ergodic phase, with the time-translational 
symmetry and the relation [d t ACij(t — t')} Cij(t — t') = d t Mij(t — t'), we can verify that the 
fluctuation-dissipation theorem indeed holds. 



From Eqs @ and ([]), a straightforward generalization of ref. ||'s scheme for the p-spin 
model leads to a closed set of integro-differential equations for the correlation functions. In 
the present case, due to the sequence-translational symmetry of ring polymers, it is more 
transparent to write the equation in the Rouse mode representation. Namely, the Fourier 
transform of the correlation function with respect to the sequence satisfies the equation 

T^dtCpit) + c p {t)/c p {0) +nf dt'd t Cp{t - t') [mo(t') - m p {t')\ = 0, (7) 

where c p (t) = FT( i _ ? -)_ >p [C'.y(£)] is the Rouse mode correlation function and m p {t) = 
FT(j_j)^p [My(t)] is a nonlinear function of all of the correlations {c p >}. The equal time 
correlation function c p (0) obeys 

Cp(O)- 1 = FT (i-jy^Kid - /i[m Q (0) - m p (0)] = k p . (8) 

If we introduce the normalized Edwards-Anderson (EA) Rouse-Zimm order parameters 
q p = lim^oo Cp(t)/c p (0), the self-consistent equation for q p becomes 

<?p/(1 - Q P ) = Ws) t m o(oo) - m p (oo)] = F p} (9) 



which is isomorphous to the equations for the MCT of structural glasses [FSj. Notice that 
if F p ^ the RZ modes have a static offset indicating the trapping in a local minimum 
and that all p modes are coupled in evaluating m p (oo). Thus, clearly, the dynamic glass 
transition, if any, should occur simultaneously for all modes in this analysis. In general, the 
dynamic glass transition can be obtained by two steps. We first solve eq.(||) for c p (0) and 
then look for a non-trivial solution of eq.(||). Assuming c p (0), a purely static quantity, does 
not exhibit any singular behavior |||18]]> we use an unperturbed confined Rouse value (i.e. 
FT (i-j)^ P [Kij}) for c p (0) here avoiding the first step. 

Fig.l shows EA parameters function of /i with p — 1, 50, 200 and 512 for a 1024- 

mer with parameters a = 1 and rj = 0.8. The EA parameters indeed exhibit a discontinuous 
(first order) transition at a critical value denoted by /xa for all p modes. Increasing the chain 
length up to 32768, we numerically show that Ta defined by /iA = (6//cbTa) 2 converges to a 



finite value Ta = 0.36(see Fig. 2), which is in good agreement with the value 0.2926 estimated 
by the static replica theory for essentially the same model [7J]. This value depends little on 
the choice of rj and thus the agreement with the replica result is quite robust; with rj = 0.6 
and rj = 1.0 we found Ta = 0.286 and Ta = 0.326, respectively. 

Scanning a wide range of chain lengths, iV, the self-consistent equations yield two different 
transitions depending on chain length and state of collapase (Fig. 2). For collapsed phase 
(e.g., 7] = 0.8), with short chain length, the EA parameters first achieve a small nonzero value 
(of order 10 -2 ) giving a weakly nonergodic phase at a parameter value /i c (dashed curve in 
Fig. 2). A stronger nonergodic phase, in which EA parameters become of order unity (/m.), 
appears for collapsed state with chain length longer than 128 (solid curve in Fig. 2). /i c 
and fi\ coalesce at iV ~ 512 above which no weakly nonergodic phase exists. Under B 
solvent conditions (i.e., 77 = 0), on the other hand, we do not find any second transition 
and fi c continues to increase with N. Since the statistical dynamical theory presented here 
is inherently appropriate for sufficiently large systems, we interpret the transition at fi^ 
to correspond to the dynamical transition found by the replica approach |J|7j . The weak 
transition (/i c ) found here only for short chains is fragile and could be an artifact of the 
model and we postpone its detailed interpretation. 

Fig. 3 plots q p = lim^oo c p (t) as a function of p as well as its inverse Fourier transformed 
frozen static displacement Qi-j. Since the p = mode is purely diffusive, we drop this 
component. Roughly, the frozen fluctuations q p and Qi-jS vary with p and i — j in a manner 
quite similar to the equal time fluctuations c p (0) and Cj-j, respectively. Although the 
frozen modes are quite localized in the bead representation, the long wavelength fluctuations 
corresponding to small p are considerable. 

Which mode is dominantly responsible for the transition? The reduction theorem of 



Gotze ||T8[ suggests that within MCT only one particular mode causes the instability of the 
nonergodic glass phase and the eigenvector of the stability matrix S pp ' = dF p / dq p / with the 
largest eigenvalue corresponds to the most dangerous mode. Dashed curves in Fig. 3 show 
the right-eigenvector of this dangerous mode in the Rouse and bead representation, v p and 



Vi-j-, have similar behavior to c p (0) and Ci-j, respectively. 

Eq.(0) can directly be integrated to get the explicit time dependence of c p {t). As in 
the MCT for structural glasses, critical slowing down is expected to be found as Ta is 



approached. Details will be given elsewhere |]TT|j . 



As mentioned above, the analysis of Roan and Shakhnovich yielded no singular behavior 
in the relaxation spectrum. There are several points in their analysis which differ from the 
present one but primarily their analysis is a lowest order perturbation with respect to bij. 
However, it is clear from our argument that the self-consistent nature of the MC calculation 
is essential for the dynamic glass transition and obviously this is included at a higher order 
in perturbation. This is a principal reason for their discordant result. Thirumalai et al.'s 
approach is quite similar in spirit to ours, but they analyze an uncoupled mode equation 
analogous to eq.(|7|) which yields a scale-dependent transition temperature. We believe that 
it is more accurate to include the coupling of different Rouse modes, forcing them to freeze 
simulutaneously. The consistency with the static replica calculation buttresses the belief. 

We emphasize again that the dynamic glass transition studied here by elementary MCT 
is, in reality, a crossover being smeared out due to entropic droplets as discussed in refs. ||[7] . 
Thus, Ta should not be viewed as a crisp phase transition, but a characteristic temperature 
at which the nature of chain dynamics changes qualitatively from renormalized free chain 
dynamics to activated escape from traps. Above Ta, the so-called cage effect due to the 
mode coupling addressed here renormalizes the Rouse relaxation through internal friction 
yielding slow dynamics analogous to the a relaxation of supercooled liquids, while below 
TA(but still above the static glass transition temperature) escape from localized free energy 
minima by (local) thermal activations, i.e. entropic droplets, controls the dynamics. 

In applying these results to protein folding we need to consider that proteins are not 
entirely random but have evolved to satisfy the minimal frustration principle |]J so as to avoid 
being trapped in nonnative local minima. The sequence is designed so that energetically 
favorable but structurally incorrect amino acid contacts are minimized. Thus, effects of 
the random heteropolymer's interaction considered here are superimposed on the flow of 
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configurations through a global funnel leading to the native structure |T|. Qualitatively, the 
folding transition temperature Tp should be larger than the static glass transition Tk for fast 
folding and may perhaps be even larger than Ta. If this is the case, folding can be viewed 
as a diffusion process in an order parameter space, where slow dynamics renormalized by 
mode coupling controls the configurational diffusion rate. 
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Figure 1: Edwards-Anderson order parameters q p for p = 1,50,200, and 512 as a function 
of fj, = (/3b) 2 . Parameters used are N = 1024, r = 1, a = 1 and n = 0.8. 
Figure 2: /xa = (b/k B TA) 2 and /x c as a function of chain length N. /xa (solid curve) and \ic 
(dashed curve) are actually computed at chain length TP with p — 5, • • • , 15. i] = for 
solvent and rj = 0.6, 0.8, and 1.0 for collapsed state. Parameters used are T = a = 1. 
Figure 3: Edwards- Anderson order parameters q p and Qi-j as a function of p and % — j 
(solid curves) and vector components v v and v^j of the dangereous mode in p and % — j 
representations (dashed curves). Parameters used are the same as Fig.l and /j is slightly 
above /xa- 
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